Self-organization in systems of self-propelled particles 
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We investigate a discrete model consisting of self-propelled particles that obey simple interaction 
rules. We show that this model can self-organize and exhibit coherent localized solutions in one- and 
in two-dimensions. In one-dimension, the self-organized solution is a localized flock of finite extent 
in which the density abruptly drops to zero at the edges. In two-dimensions, we focus on the vortex 
solution in which the particles rotate around a common center and show that this solution can be 
obtained from random initial conditions, even in the absence of a confining boundary. Furthermore, 
we develop a continuum version of our discrete model and demonstrate that the agreement between 
the discrete and the continuum model is excellent. 



Self-organization and pattern formation in systems of 
self-propelled entities are ubiquitous in nature. Exam- 
ples can be found in a variety of fields and include ani- 
mal aggregation J|] , traffic patterns || and cell migration 
j|. Recently, the problem of flocking, in which a large 
number of moving particles (e.g. fish or birds) remain co- 
herent over long times and distances, has attracted con- 
siderable attention. A simulation of a simple numerical 
model by Vicsek et al. [0 in which each particle has a 
constant identical speed and a direction of motion that 
is determined by the average direction of its neighbors, 
revealed that an ordered phase exists, even in the present 
of noise and disorder. A subsequent analysis of a contin- 
uum model by Toner and Tu || investigated this ordered 
phase further and derived conditions for its stability. 

These models have in common that the flocks have 
infinite extent and, in simulations, fill the entire compu- 
tational box. In reality, however, flocks have a finite size, 
with its density dropping sharply at the edge of the flock 
||. In this Letter we present a discrete model consist- 
ing of self-propelled interacting particles that obey simple 
rules. We show that self-organization in our model leads 
to coherent localized states in one dimension (ID) and 
in two dimensions (2D) that are stable in the presence 
of noise and disorder. Furthermore, we present a con- 
tinuum version of our discrete model which is obtained 
by coarse-grain averaging the discrete equations. The 
continuum flock solutions in ID agree very well with the 
discrete solutions and are characterized by having a finite 
extent with densities that drop off sharply at the edges. 
In 2D, we focus on a vortex state in which the particles 
are rotating around a common center and show that the 
discrete model and the continuum model agree well. 

Our discrete particle model consists of N particles with 
mass mi, position Xi and velocity Uj. Each particle ex- 
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FIG. 1. A coherent moving flock in the one-dimensional 
version of the model with parameters (all with arbitrary units) 
m=l, a=0.5, /3=1, Ca=0.45, L=60, CV=2 and l r =20. The 
solid circles correspond to the solution of the discrete model 
for 7V=200 and every 10th particle displayed. The solid line 
shows the solution of the continuum model. 

periences a self-propelling force ft with fixed magnitude 
a. To prevent the particles from reaching large speeds, 
a friction force with coefficient j3 is introduced. In addi- 
tion, each particle is subject to an attractive force which 
is characterized by an interaction range l a . This force is 
responsible for the aggregation of the particles and corre- 
sponds in animal aggregates to an awareness of the posi- 
tion of surrounding animals. To prevent a collapse of the 
aggregate, a shorter-range repulsive force is introduced 
with interaction range l r . Thus, the governing equations 
for each particle is 

rriidtVi = af t - fiv t - VU (1) 
d t Xi = Vi (2) 
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FIG. 2. The size L of a flock as a function of N with (solid 
line) and without (dashed line) an additional short-range 
hard-core potential. The inset shows the density of flocks in 
the presence of the short-range hard-core potential for differ- 
ent N. Parameter values: m=l, a=0.5, ,9=1, C a =0.6, Z a =40, 
CV=2, l r =20, C hc =l and l hc =W. 



We have checked that our qualitative results are indepen- 
dent of the explicit form of the interaction potential and 
we have chosen here an exponentially decaying interac- 
tion: 



U = exp(-|xj ~Xj\/l a ) 



C r exp(— \xi — Xj | / l r ) 



(3) 



where C a , C r determine the strength of the attractive 
and repulsive force respectively. The direction of the self- 
propelling force can be chosen along the instantaneous 
velocity vector or, similar to the numerical model of Ref. 
[^J, can be determined by aligning it with the average 
velocity direction of the neighboring particles: 



fi = Hi without averaging 
Vj exp(— \xi — Xj\/l c ) with averaging (4) 



(5) 



where l c is a correlation length. 

Let us now present our simulation results of the dis- 
crete model. The model was integrated by solving Eqs. 
|l|,H using a simple Euler integration routine with timestep 
At = 0.2. The simplest coherent localized solution in our 
model is a ID flock in which all particles move with con- 
stant velocity v = a/ (3. An example of this solution is 
presented in Fig. |l| where we have plotted, as solid cir- 
cles, the density defined as pi = 2/(x(i + 1) — x(i — 1)) 



as a function of the position of the particle. The density 
can be seen to drop abruptly to zero at the edge of the 
flock. We have checked that this solution is stable in the 
presence of moderate amounts of noise (added to Eq. |^) 
and of disorder in the parameters. 

Further simulations revealed that increasing the num- 
ber of particles does not change the shape of the density 
function and that the total size of the flock reaches a con- 
stant value. This is illustrated in Fig. [2] where the dashed 
line represents the size of the flock as a function of the 
number of particles. This obviously unrealistic behavior 
of the model is due to the soft-core repulsive force which 
allows the particles to be very close. Our model can easily 
be extended to incorporate hard-core repulsive forces. In 
fact, a flock with a size that scales approximately linearly 
with the total number of particles can be obtained by in- 
troducing a hard-core repulsive force in Eq. [l] 0. The 
specific form of the hard-core potential is not important 
and we have chosen 



I he) 5 



j I — ^hc 



(6) 



U hc = | 

In Fig. H we show, as a solid line, the size of the flock in 
the presence of this additional repulsive force as a func- 
tion of N while in the inset we show the corresponding 
density of the flock for different N. Naturally, the force 
has only an effect when the interparticle separations are 
smaller than lh c - Hence, for small N the flock solution 
is unaffected by the additional force. As N is increased 
and the interparticle spacing becomes smaller than Z^ c , 
the particles in the center of the flock are pushed apart. 
For large N, the flock reaches a constant density in its 
center and its size scales linearly with N. 

Let us now turn to 2D where we have obtained several 
flocking states. One, not shown here, is the equivalent of 
our ID flock: all particles are moving in the same direc- 
tion with \vl\ = a/ [3. The particles arrange themselves on 
a disk and this aggregate is stable under small amounts 
of noise and disorder. The solution which we will focus 
on here consists of a vortex state where the particles ro- 
tate around a common center and which is common in 
fish schools ||, bacterial colonies and amoeba aggre- 
gates [jlO^I . This solution has been observed previously in 
models of self-propelling particles but only in the pres- 
ence of a confining boundary Tl|Jl2] ] or when a rotational 
chemotaxis term is invoked K ] . 

In our model, the vortex solution can be obtained from 
a wide variety of initial conditions including one in which 
all particles are randomly placed on a disk with speed 
a/ (3 and random initial velocities. A typical evolution of 
the particles starting from this initial condition is shown 
m Fig. | This fi gure was obtained in the absence of 
the velocity averaging term and illustrates that in this 
case some particles move clockwise while the others ro- 
tate counter clockwise. When the velocity averaging term 
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FIG. 3. Snapshots of A?=200 particles for the parameter 
values a=10, f3=l, l c = 0, C a =0.4, Z a =40, C r =l and l r =20. 
As initial condition, the particles are placed at random on 
a disk with velocities that are constant in magnitude (a/ (3) 
but random in direction (a). After an initial transient (b, 20 
iterations and c, 50 iterations), a stable rotating vortex state 
is formed (d, 300 iterations). The bar indicates the attraction 
length l a - The position of each particle is denoted by a solid 
circle and the velocity as a line starting at the particle and 
pointing in the direction of the velocity. 

is included the final vortex state consists of all particles 
turning the same way. An example of this case is shown 
in the inset of Fig. ^. In both cases, the speed of each 
particle was found to be sharply peaked around \ v\ — a/ f3 

The average size of the vortex remains constant in time 
as shown in Fig. ^ where we have plotted the average 
density (obtained by averaging over 10 6 iterations) of the 
vortex structure. Fig. ^ displays several remarkable fea- 
tures. First, there is a well-defined core which remains 
void even for extended simulation runs. Using different 
parameter values however, one can also produce a vortex 
without a core. Second, as in the one dimensional flock, 
the density does not decay smoothly to zero at the edges. 
Instead, it increases at both the inner and outer edge 
of the aggregate and then drops abruptly. Qualitatively 
similar vortex solutions were found when an additional 
hard-core repulsion like the one discussed above is added. 

As in the case of traffic models jl3| , it is useful to de- 
velop a continuum version of our model. To this end, we 
simply coarse-grain average the equations which results 
in, after dividing by a common factor of p, 

d t v + (v- V)u = af-j3v + G (7) 

together with the conservation equation for the density 
P 
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FIG. 4. Average density of a rotating vortex state in 
the discrete model (solid symbols) and the continuum model 
for the parameter values A r =400, m=l, a=10, /3=1, l c = 0, 
C a ~0.5, i a =30, C r =l and 1, -20. The inset shows a snapshot 
of the discrete model simulation with the bar corresponding 
to la- As initial condition we used a vortex obtained with 
l c — 4 which ensured that the angular velocity of all particles 
has the same sign. 

d t p + V • (vp) = 0. 
The interaction force is given by 

aw- f ro 

and the self-propulsive force direction is given by either 
f{x) = J p{x')v(x')exp{-\x-x'\/l d )dx' (9) 

in the velocity averaging case or simply by / = v/\v\ 
otherwise. 

A comparison between the discrete model and the con- 
tinuum model can be carried out for the solutions pre- 
sented here. For a ID flock, we have simply fi = 1, and 
the solution in the continuum model is given by G = 
or 

JfWU^^U.D (10) 

where D is a constant determining the total number of 
particles. Since the sought after solution has a finite ex- 
tent with a discontinuity at the edge where the density 
drops to zero we discretize the integral using M points 
and discretization step Ax. The last point corresponds 
to the edge of the flock. The resulting linear set of M 
equations for p is easily solved using standard linear al- 
gebra packages, and Ax was varied until the slope at the 
center of the flock vanished. The result, with D chosen 
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FIG. 5. Two different solutions of the continuum model 
for the parameters a//3=10, C a =0.7, £ a =40, C r =2, 1,-20 and 
7V=200. 

such that J p(x)dx — N , is shown in Fig. [I] as a solid 
line. The density profile in the continuum model is dis- 
continuous at the edge and agrees well with the profile 
obtained in the discrete model. Note that since the equa- 
tions are linear in p it is not surprising to find that the 
density profile of the discrete model did not change as 
the number of particles is increased. Clearlythe simple 
coarse-grained averaging procedure is not adequate for 
the hard-core potential case, where higher order terms in 
the density are important. 

The continuum model can also be used to find the vor- 
tex solution. To this end, we use the fact that all particles 
undergo approximately a circular motion with constant 
speed {a/ [3). Thus, a continuum vortex solution can be 
found by requiring that the force G is centripetal: 

/>27T />0O 

/ d<P drp(r)U(r,<P) = D - (a / f3) 2 \n(r) (11) 
Jo Jo 

After performing the integration over <f>, the remaining 
integral was discretized as in the ID case. The resulting 
matrix was solved for p(r) and used in a Newton solver 
that searched for the size of the hole and the overall size 
of the vortex (i.e. discretization Ax) with a condition 
for smoothness of the solution at both discontuous edges. 
This condition simply consisted of requiring that the first 
and last point can be obtained by linear interpolation 
from its two neighboring points. We have checked that 
the solutions we obtained are converged by increasing the 
number of discretization points from 80 to 1480. In Fig. 
^ we compare the discrete solution to the one obtained 
by Eq. [H] where the integration constant D was varied 
until J p(x)dx — N. Again, the continuum profile is dis- 
continuous at the edges and the agreement between the 
continuum profile and the discrete profile is very good. 

The continuum equation can be used to explore the 
(large) parameter space more efficiently. An example of 
such an exploration is shown in Fig. ^| where we plot 
two different solutions found by our Newton solver for 



the same model parameters and total number of parti- 
cles (./V=200) but with different integration constant D. 
Preliminary simulations of the discrete model show that 
the solution with the larger core is unstable and a for- 
mal stability analysis of the continuum solution will be 
carried out in the future. 

In this Letter we have presented a model for localized 
aggregates and flocks. Our model contains very simple 
rules and can be straightforwardly extended to incorpo- 
rate additional and different types of interactions. For 
example, the forces that maintain bacterial aggregates 
are believed to be short-range adhesion forces together 
with a short-range hard-core repulsion. The investiga- 
tion of these types of interactions will be the subject of 
future work. Finally, it would be interesting to com- 
pare our results to animal flocks. Unfortunately, such 
a comparison is currently not possible since not enough 
quantitative data is available. 
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